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Abstract 

We investigate the electronic stress tensor of the hydrogen molecular ion H^" for the ground state 
using the exact wave function and wave functions approximated by gaussian function basis set 
expansion. The spatial distribution of the largest eigenvalue, corresponding eigenvectors, tension 
and kinetic energy density are compared. We find that the cc-pV6Z basis set gives the spindle 
structure very close to the one calculated from the exact wave function. Similarly, energy density 
at the Lagrange point is very well approximated by the cc-pV5Z or cc-pV6Z basis sets. 
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I. INTRODUCTION 



There have been many studies of stress tensors in quantum systems as is found in Refs [T]- 
, and even earlier f2U\ |2"T] . Several ways to define and apply the stress tensor in quantum 
systems are proposed in the literature. In this paper, we consider the electronic stress tensor 
and its application to chemical bonds and reactions as has been studied in Refs. pHl EEl 
IH [T6l [TH [22H21]. In these works, electronic stress tensors which measure effects caused by 
internal forces acting on electrons in molecules play a central role, in contrast to Ref. [I] and 
followers who focus on the stress tensor which is associated with forces on nuclei. 

Recently, the authors studied the electronic stress tensor of hydrogen molecular ion 
using its exact wave function [25] (Refs. [26l E7j investigated stress tensor field of a Hj" 
molecule but with different definitions from ours). By using the exact wave functions, we 
were able to investigate the stress tensor and other associated quantities in a very accurate 
manner. However, since we can only obtain approximate wave functions for almost every 
molecular system, it is very important to study how the stress tensor of deviates from 
the exact form as we change the level of approximation of the quantum chemistry calculation 
to derive wave functions. This is the purpose of this paper. We believe this study helps 
to solidify the basis of the stress tensor analysis of molecular systems and open up areas of 
further application. 

This paper is organized as follows. In Sec. [TT| we introduce quantities which are used in 
our stress tensor analysis. We also explain the exact and approximate wave functions we 



use. In Sec. Ill, we compare the quantities obtained using the exact wave functions and 



approximate wave functions. Sec. IV gives our conclusion. 



II. CALCULATION METHODS 

A. Electronic stress tensor, tension and energy density 

We first show our expression of the electronic stress tensor t$ constructed from the 
wave function ip. That is, 



2 



where {i,j} = {1, 2, 3} denote spatial coordinates, m is the electron mass and c.c. stands for 
complex conjugate. The expression for many-electron systems, that is for almost every other 
molecules, is found in Ref. [10]. There, the stress tensor has been derived by field theoretic 
method which is applicable to many-particle systems. Here, for illustrative purpose, we 
derive it using the Schrodinger equation for a single electron in H^~. We adopt the adiabatic 
approximation in this paper. 

The time-dependent Schrodinger equation reads 

where the potential energy V is assumed to be real. It is well known that the continuity 
equation dn/dt + V • j = holds by defining probability density n = \ip\ 2 and probability 
flux j = (7&/m)Im(^>*V'0)- Now, the equilibrium equation for the electronic stress tensor is 
obtained from an equation of motion of j. Namely, using Eq. (|2j) , time derivative of j can 
be written as 

% = (V^ • VV - V*VV 2 V + c.c.) - Kw. 
ot 4m z ' m 

For the steady state which is the case of molecular systems we consider, since dj/dt = 0, 
we obtain the equilibrium equation 

djrfj + F L>i = 0, (3) 

where = —\ip\ 2 diV is the Lorentz force generated by nuclei. For later convenience, we 
here define the first term in the equation as tension F T ^ = djrfj. From this equation, we see 
that t$ is the stress caused by purely quantum mechanical effects. The tension is written 
explicitly as 

h 2 

F T = (VV> • VV - V^*VV 2 ^ + c.c.) . (4) 

4m v ' 

Since this force, which cancels the classical Lorentz force at each point in space (Eq. d3])), 
also expresses purely quantum mechanical effects, it is considered to carry some information 
of a chemical bond. Actually, in Ref. [16] , the point where the tension vanishes is defined as 
"Lagrange point" and it is proposed that we may calculate energy density [defined below, 
Eq. (|6f ] at that point to define bond order. This new definition turns out to carry very nice 
features of bond order [16l 



We now give our definition of energy density e T constructed from the stress tensor. In 
Ref. [10J, it has been proposed to define energy density from the trace of the stress tensor 
as 

i 

h 2 

= — U*v 2 ij - v^* ■ Vi) + c.c.) . (6) 

8m 

Note that this definition gives correct total energy when integrated over the whole space 
and the virial theorem is applied. 

We finally describe the kinetic energy density ny, which is given by 

n T = -^ l -{rV^ + c.c). (7) 

Using kinetic energy density, we can divide the whole space into the region with > 
0, where the electron classical motion is allowed, ut < 0, where the classical motion is 
forbidden, and ut = 0, the boundary between them. These regions are respectively called the 
electronic drop region Rp, the electronic atmosphere region Ra and the electronic interface 
S [ID]- The S can give the effective shape of the molecule and is therefore an important 
region in particular. 

The numerical calculation of these quantities are performed by the MRDFT code devel- 
oped by our group [28J. 



B. Exact and approximate wave functions 



We calculate quantities defined in Sec. II A using the exact wave functions and approxi- 



mate wave functions expanded in Gaussian basis set. We use the exact wave function which 
is computed in Ref. [29]. The calculations of the stress tensor and tension using this wave 
function have been reported by the authors in Ref. [25J . In the present paper, we focus on the 
ground state and fix the internuclear distance to be the equilibrium distance R e = 2.0bohrs. 
(To be more precise, R e = 1.9972 bohrs [30] but it does not make practical difference in our 
argument below.) 

For the calculation of approximate wave functions, we use Gaussian 09 code [21] to 
perform the Hartree-Fock calculation with the following basis sets: STO-3G [32J, 6-31G 
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[33] . 6-311G [31] and cc-pVxZ (x = D, T, Q, 5 and 6) [351 EE]- For 6-31G, we also consider 
adding polarization and/or diffuse functions (6-31++G, 6-31G** and 6-31++G**) [33 EH], 
and similar for 6-311G (6-311++G, 6-311G** and 6-311++G**) [Ml ES]- The cc-pVxZ 
basis sets include polarization functions in their definition and we consider adding diffuse 
functions to each of them (aug-cc-pVxZ) [361 EE]. We use Cartesian functions (i.e. 6 d, 10 
f, 15 g and 21 h functions) for all of the basis sets. Namely, we consider 19 types of basis 
set. We note that since is a one-electron system, the accuracy of the calculation depends 
only on the basis set. 

III. RESULTS AND DISCUSSION 

In this section, we calculate the quantities introduced in the previous section for the 
exact and approximate wave functions of Hj. We compare among them and investigate 
how the results differ from those of the exact wave functions depending on the basis sets. 
In the following figures, we take the origin of the coordinate to be the midpoint of the two 
H nuclei and z-axis to be the internuclear axis. We locate H atoms at (z,x) = (—1.0,0.0) 
and (1.0,0.0). Although this system has the rotational symmetry around the z-axis and 
reflection symmetry for the x = plane, we sometimes draw figure for the regions with 
x < or z < to facilitate the intuitive understanding of chemical bond of H^~. 

A. Electron density 

Before we turn to discuss the electronic stress tensor, tension and so forth, it may be 
useful to see how the electron density varies with respect to the basis set we use. For this 
purpose, we plot the spatial distribution of relative error of the electron density An in Fig. [Tj 
This is computed as 

An = na PP~ Ue \ (8) 

where n app and n ex are the electron density from approximate wave functions and the exact 
wave function respectively. When we use cc-pV5Z or larger, the electron density is calculated 
with less than 1% relative error in the region shown in the figure (\x\ < 3 and \z\ < 3). The 
best performance in this region is achieved by cc-pV6Z, which has 0.13% error, rather than 
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aug-cc-pV6Z, which has 0.44% error. We can also confirm that adding polarization functions 
makes wave functions closer to the exact wave function especially in the region between nuclei 
whereas adding diffuse functions make little difference. 

B. Eigenvalue and eigenvector of the stress tensor 

We here discuss the electronic stress tensor, defined by Eq. ([I]). As is often done for the 
stress tensor in general, we examine its largest eigenvalue and corresponding eigenvector. 
The sign of the largest eigenvalue tells whether electrons at a certain point in space feel tensile 
force (positive eigenvalue) or compressive force (negative eigenvalue) and the eigenvector tells 
direction of the force. They are shown in Fig. [2] for the exact wave function and approximate 
wave functions. In Fig. [3j we also plot the relative error of the largest eigenvalue Ae 

Ae = 6(W ~ eex , (9) 

where e app and e ex are the largest eigenvalues from approximate wave functions and the 
exact wave function respectively. To show the details of the direction of the eigenvector v eig , 
we plot 

e eig = arctan f J , (10) 

in Fig. [4j This quantity ranges — ir/2 < 9 eig < n/2 and, in our z-x plane, positive slope 
takes positive value (shown in red) and negative slope takes negative value (shown in blue). 

Let us first point out the properties found in the case of the exact wave function [25] . The 
region of the positive eigenvalue occupies most of the space between the nuclei. In detail, 
the positive eigenvalue region is bounded by a closed sphere-like surface that touches two H 
nuclei. Also, in that region, the eigenvectors forms a bundle of flow lines that connects the H 
nuclei. This region, called "spindle structure" [TB], is clearly visible using the analysis with 
the exact wave function. Such structure expresses the correct directionality of the chemical 
bond and the positive eigenvalue of the stress tensor, implying a tensile stress, is considered 
to well characterize the covalent bond. 

As for the approximate wave functions, we see the region with positive eigenvalue (shown 
in red) spreads between the nuclei and also flows of the eigenvectors connecting the H nuclei. 
We can say that qualitative features, especially the spindle structure, are reproduced by the 
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approximate wave functions. In detail, however, there are some differences. For one thing, 
the positive eigenvalue region in approximate cases surrounds the nuclei and it does not 
touch them (it is not simply connected in contrast to the case of the exact wave function). 
To put it another way, the positive region spreads outside the internuclear region and the 
negative region exists within the internuclear region. For the case of the exact wave function, 
there is only one zero surface of the eigenvalue and the whole space is divided into just a 
positive region and a negative region, but we see multiple zero surfaces for the case of the 
approximate wave functions in particular nearby and away from the nuclei. 

Comparing among the approximate wave functions, we see that the cc-pV6Z seems to 
quite closely reproduce the exact case. The cc-pV5Z also reproduces well the exact case. 
We may say that we need to approximate the wave function with less than 1% error to have 
a good stress tensor profile. Similarly to the case of the electron density, the addition of 
the polarization functions improves the shape of the spindle structure whereas the diffuse 
functions do not. This is reasonable because the spindle structure is closely connected to 
the nature of bonding. 

It may be interesting to point out that the pattern of eigenvectors nearby the nuclei (that 
they point radially from the nuclei) does not show much difference between the exact and 
approximate wave functions even for the basis sets which are not so large. Furthermore, such 
resemblance holds even at places where the signs of the eigenvalues are different between 
the exact and approximate cases. This is somewhat surprising feature and may be just a 
coincidence but we do not have an explanation. 



C. Tension 

Next, we examine the tension field computed as Eq. Q. This is in Fig. [5] where directions 
are shown by arrows and norm is shown by a color map. For the details of the direction, we 
plot 

Oten = atan2 (F T)X , F T>Z ) , (11) 
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in Fig. [6j 1 This is the angle between the positive z-axis and the point (F T)Z , F T>X ), with the 
range — tt < 9 ten < n. In our z-x plane, the angle is positive for counter-clockwise angles 
(upper half-plane, x > 0), and negative for clockwise angles (lower half-plane, x < 0). 

When we compare the tension of the exact and approximate wave functions, it seems that 
6-311G and cc-pVa;Z (x = D, T, Q, 5 and 6) do a good job in reproducing the exact result. 
In other words, compared with the case of the eigenvalue of the stress tensor, the tension 
can be computed rather accurately with smaller basis sets. This reflects the fact that the 



tension is just the opposite of the Coulomb force F L = —\ip\ 2 'VV (see Sec. II A) and does 
not have much information as the stress tensor. 

The most important feature of the tension is the Lagrange points, where the tension 



vanishes, as mentioned in Sec. II A In the case of Hj , the Lagrange point turns out to be 
the midpoint of two H nuclei, (0.0, 0.0). Actually, we see the region with small value of \F T \ 
(expressed in white) in the neighborhood of the origin for the exact wave function. Visual 
inspection of the figure tells us that such region is found in 6-311G, cc-pVTZ, cc-pV5Z and 
cc-pV6Z but not for cc-pVDZ and cc-pVQZ. Hence, relation between the area of the region 
with small \F T \ and the level of the basis set seems to be somewhat irregular. In the case 
of Hj, we can find a zero point in every basis set if we look at the data closely, but it may 
happen in general that the Lagrange points are found for smaller basis sets whereas not for 
larger basis sets. 

For the direction of the tension vector, as is found in the case of the eigenvector of the 
stress tensor discussed above, it is interesting that approximate wave functions give similar 
pattern to that of the exact wave function near the positions of the nuclei. 



1 This function can be written using the standard arctan function as follows: 

arctan(a;/z) z > 

7r + arctan(a;/z) x > 0, z < 

atan2(x, z) = ( —ir + arctan(a;/z) x < 0, z < (12) 

tt/2 x>0,z = 

-7I-/2 x<0 1 z = 
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D. Kinetic energy density 



We now examine the kinetic energy density defined as Eq. (J7|). This is plotted in Fig. [7j 
The zero surface of the kinetic energy density can be considered as an effective surface of a 
molecule and, for the exact wave function of the Hj" molecule, it is a closed surface which 
includes the H nuclei. This feature is reproduced by every approximate wave function but 
if we want to reproduce the detailed pattern of the curvature of the surface, it seems that 
we need cc-pV5Z or cc-pV6Z. 

In contrast to the case of the zero surface of the largest eigenvalue of the stress tensor 
discussed above, there is no extra zero surface around the nuclei for most of the basis sets 
(exceptions are 6-31G** and 6-31++G**). 



E. Quantities at Lagrange point 



As is mentioned in Sec. II A the Lagrange points are proposed to be useful to charac- 
terize chemical bond [TB]. In particular, new bond order has been defined as the energy 
density (Eq. (J6|) at the Lagrange points [TBI E21 E3]- Therefore, it is important to check 
the energy density at the Lagrange point is reproduced by approximate wave functions. We 
also compare the largest eigenvalue of the stress tensor. 

In our case of the molecule, we have a Lagrange point at the midpoint of two H nuclei. 
We calculate the electron density, the largest eigenvalue of the stress tensor and energy 
density at the Lagrange point for the exact and approximate wave functions (Table [T]) and 
also relative errors between those quantities derived from the exact and approximate wave 
functions (Table [TT] and Fig. [8]). Roughly speaking, the largest eigenvalue and energy density 
at the Lagrange point are calculated with about 10 times larger relative error than that of 
electon density. We see that cc-pV5Z and cc-pV6Z can reproduce almost same results (less 
than 1% relative error) as the exact wave function. 



IV. CONCLUSION 



In this paper, we investigated the electronic stress tensor of the hydrogen molecular ion 
for the ground state using the exact wave function and wave functions approximated by 
gaussian function basis set expansion. We compared the spatial distribution of the largest 

9 



eigenvalue of the stress tensor, corresponding eigenvectors, tension and kinetic energy density 
between the exact wave functions and approximate wave functions using 19 types of basis 
sets. We also compared the energy density and the largest eigenvalue at the Lagrange point. 

We found that the approximate wave function using the cc-pV6Z basis set gives the 
spindle structure, which is the key structure expressing the covalent bond, very close to the 
one calculated from the exact wave function. We also found that the energy density at the 
Lagrange point, which is the quantity used to define bond order, calculated by the exact 
wave function is very well approximated by the calculation using the cc-pV5Z or cc-pV6Z 
basis sets. 

It is a good news that the electronic stress tensor and associated quantities can be com- 
puted accurately with commonly used basis sets which are stored in Gaussian program. Next 
step would be to compute the stress tensor for other molecules accurately and investigate 
their spindle structure in detail. 
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TABLE I: Electron density, the largest eigenvalue of the stress tensor and energy density at the 
Lagrange point (z,x) = (0.0,0.0). All the values are multiplied by 10 2 . 



Basis set 


Electron density 


Largest Eigenvalue 


Energy density 


STO-3G 


6.801 


5.943 


-1.063 


6-3 1G 


6.916 


6.219 


-1.246 


6-31++G 


6.960 


6.128 


-1.271 


6-31G** 


9.430 


5.947 


-3.674 


6-31++G** 


9.463 


5.838 


-3.694 


6-311G 


7.398 


4.233 


-2.231 


6-311++G 


7.407 


4.228 


-2.238 


6-311G** 


9.725 


3.287 


-4.385 


6-311++G** 


9.720 


3.292 


-4.381 


cc-pVDZ 


9.685 


4.671 


-3.821 


aug-cc-pVDZ 


9.562 


4.552 


-3.688 


cc-pVTZ 


9.948 


3.356 


-4.756 


aug-cc-pVTZ 


9.943 


3.318 


-4.768 


cc-pVQZ 


9.880 


4.426 


-4.222 


aug-cc-pVQZ 


9.863 


4.477 


-4.156 


cc-pV5Z 


9.902 


3.995 


-4.445 


aug-cc-pV5Z 


9.898 


4.018 


-4.416 


cc-pV6Z 


9.902 


4.028 


-4.433 


aug-cc-pV6Z 


9.901 


4.037 


-4.423 


Exact 


9.903 


4.019 


-4.443 
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TABLE II: The relative error of the electron density, the largest eigenvalue of the stress tensor and 
energy density with respect to the exact wave function at the Lagrange point (z,x) = (0.0,0.0). 
The bar plot for this table is shown in Fig. [8j 



Basis set 


Electron density 


Largest Eigenvalue 


Energy density 


STO-3G 


-3.132x10' 


-i 


4.787x10" 


l 


-7.607x10" 


-l 


6-3 1G 


-3.016x10" 


-i 


5.473x10" 


-i 


-7.196x10" 


-l 


6-31++G 


-2.972x10" 


-i 


5.247x10" 


-i 


-7.139x10" 


-l 


6-31G** 


-4.776x10" 


-2 


4.797x10" 


-i 


-1.732x10" 


-l 


6-31++G** 


-4.441x10" 


-2 


4.526x10" 


-i 


-1.687x10" 


-l 


6-311G 


-2.530x10" 


-1 


5.325x10" 


-2 


-4.980x10" 


-l 


6-311++G 


-2.520x10" 


-1 


5.185x10" 


-2 


-4.964x10" 


-l 


6-311G** 


-1.800x10" 


-2 


-1.822x10" 


-1 


-1.300x10" 


-2 


6-311++G** 


-1.853x10" 


-2 


-1.810x10" 


-1 


-1.391x10" 


-2 


cc-pVDZ 


-2.200x10" 


-2 


1.622x10" 


-1 


-1.400x10" 


-1 


aug-cc-pVDZ 


-3.448x10" 


-2 


1.325x10" 


-1 


-1.699x10" 


-1 


cc-pVTZ 


4.537x10" 


3 


-1.650x10" 


-1 


7.038x10" 


2 


aug-cc-pVTZ 


4.004x10" 


3 


-1.744x10" 


-1 


7.305x10" 


2 


cc-pVQZ 


-2.323x10" 


-3 


1.012x10" 


1 


-4.978x10" 


-2 


aug-cc-pVQZ 


-4.079x10" 


-3 


1.138x10" 


1 


-6.461x10" 


-2 


cc-pV5Z 


-9.602x10" 


-5 


-5.961x10" 


-3 


3.359x10" 


-4 


aug-cc-pV5Z 


-5.210x10" 


-4 


-3.635x10 


-4 


-6.084x10" 


-3 
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FIG. 1: The spatial distribution of the relative error between the electron density calculated from 
the exact and approximate wave functions of H^" molecule An (Eq. Q) is plotted in the plane 
including two H nuclei (located at (z,x) = (—1.0,0.0) and (1.0, 0.0) ) for various basis sets. Note 
that the scale is ten times smaller after the cc-pVTZ basis set. 
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FIG. 2: The spatial distribution of the largest eigenvalue and corresponding eigenvectors of the 
stress tensor of molecule is plotted in the plane including two H nuclei (located at (z, x) = 
(—1.0, 0.0) and (1.0, 0.0) ) for various basis sets. The black solid line shows the zero surface of the 
eigenvalue. 
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FIG. 3: The spatial distribution of the relative error between the largest eigenvalue of the electronic 
stress tensor calculated from the exact and approximate wave functions of H^" molecule Ae (Eq. Q) 
is plotted similarly to Fig. [Tj 
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FIG. 4: The spatial distribution of inclination angle 9 e i g (Eq. (|10[)) of eigenvectors measured from 
the z-axis is plotted for ranges x > and z > 0. 9 e i g is shown in units of it radian and its range is 

-7T/2 < 6 ei g < 7T/2. 
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FIG. 5: The spatial distributions of the tension F T of Hj" molecule is plotted similarly to Fig. |2| 
Normalized tension vectors F T /\F T \ are shown by arrows and the norm \F T \ is depicted by a color 
map. 



19 




0.5 1 1.5 2 2.5 3 

z (bohrs) 



1 




3 r 






2.5 


0.5 










2 - 





o 


1.5 


-a 






1 - 


-0.5 










0.5 - 


-1 




' 



1 




3 






2.5 


0.5 










2 





o 


1.5 








1 


-0.5 










0.5 


-1 








0.5 1 1.5 2 2.5 3 

z (bohrs) 




1 




3 






2.5 


0.5 










2 





o 


1.5 








1 


-0.5 










0.5 


-1 








0.5 1 1.5 2 2.5 3 

z (bohrs) 




0.5 1 1.5 2 2.5 3 

z (bohrs) 




1 




3 r 






2.5 


0.5 










2 








1.5 












1 - 


-0.5 










0.5 


-1 




o ' 



0.5 1 1.5 2 2.5 3 
z (bohrs) 



1 

0.5 



-0.5 
-1 



0.5 1 1.5 2 2.5 3 

z (bohrs) 




1 




3 






2.5 


0.5 










2 





o 


1.5 












1 


4).5 










0.5 


-1 








0.5 1 1.5 2 2.5 3 
z (bohrs) 




0.5 1 1.5 2 2.5 3 

z (bohrs) 



cc-pVDZ 



AUG-cc-pVDZ 




0.5 1 1.5 2 2.5 3 

z (bohrs) 



1 




3 






2.5 


0.5 








"e 


2 





a 


1.5 


-a 






1 


-0.5 










0.5 


-1 









1 

0.5 



-0.5 
-1 



0.5 1 1.5 2 2.5 3 

z (bohrs) 



3 
2.5 

2 
1.5 

1 

0.5 





1 




3 






2.5 


0.5 








n 


2 





o 


1.5 








1 


-0.5 










0.5 


-1 








0.5 1 1.5 2 2.5 3 
z (bohrs) 




0.5 1 1.5 2 2.5 3 

z (bohrs) 



1 

0.5 



-0.5 
-1 



cc-pVTZ 



AUG^c-pVTZ 



cc-pVQZ 



AUG-cc-pVQZ 



3 






1 




3 






i 


2.5 
2 




0.5 




2.5 
2 






0.5 


1.5 









boh 


1.5 









1 

0.5 




: 


: 




1 

0.5 






-0.5 


















U-, 



0.5 1 1.5 2 2.5 3 
z (bohrs) 



0.5 1 1.5 2 2.5 3 

z (bohrs) 




1 




3 






2.5 


0.5 










2 





o 


1.5 








1 


JX5 










0.5 


-1 




o ' 



0.5 1 1.5 2 2.5 3 
z (bohrs) 



1 

0.5 



-0.5 
-1 



0.5 1 1.5 2 2.5 3 

z (bohrs) 



cc-pV5Z 



AUG-cc-pV5Z 



cc-pV6Z 



AUG-cc-pV6Z 




1 




3 








2.5 




0.5 












2 







o 


1.5 












1 




-0.5 












0.5 




-1 










0.5 1 1.5 2 2.5 3 

z (bohrs) 



1 




3 






2.5 


0.5 










2 





o 


1.5 








1 


-0.5 










0.5 


-1 








0.5 1 1.5 2 2.5 3 

z (bohrs) 




1 




3 






2.5 


0.5 










2 





o 


1.5 








1 


4X5 










0.5 


-1 








0.5 1 1.5 2 2.5 3 
z (bohrs) 



1 

0.5 



-0.5 
-1 



0.5 1 1.5 2 2.5 3 

z (bohrs) 



FIG. 6: The spatial distribution of direction angle 6 ten (Eq. (11)) of tension vectors measured from 
the z-axis is plotted for ranges x > and z > 0. 9 te n is shown in units of ir radian and its range is 

-7T < 9 ten < IT. 
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FIG. 7: The spatial distribution of kinetic energy density is plotted similarly to Fig. [2} The black 
solid line shows the zero surface of the kinetic energy density. 
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FIG. 8: The relative error of the electron density, the largest eigenvalue of the stress tensor and 
energy density with respect to the exact wave function at the Lagrange point (z, x) = (0.0,0.0). 
See Table HH for the detailed numbers. 
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